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Abstract 

What is the response of the climate system to anthropogenic forcings? 
This question is addressed typically using projections from climate models. 
The uncertainty surrounding current climate projections has important pol- 
icy implications. Characterizing and, if possible, reducing this uncertainty 
is an area of ongoing research. We consider the problem of making pro- 
jections of the North Atlantic meridional overturning circulation (AMOC). 
AMOC projections are of interest because AMOC changes may considerably 
impact natural and human systems. Uncertainties about climate model 
parameters play a key role in uncertainties in AMOC projections. When 
the observational data and the climate model output are high-dimensional 
spatial data sets, the data are typically aggregated due to computational 
constraints. The effects of aggregation are unclear because statistically rig- 
orous approaches for model parameter inference have been infeasible for 
high-resolution data. Here we develop a flexible and computationally effi- 
cient approach using principal components and basis expansions to study the 
effect of spatial data aggregation on parametric and projection uncertainties. 
Our Bayesian reduced-dimensional calibration approach allows us to study 
the effect of complicated error structures and data-model discrepancies on 



23 



our ability to learn about climate model parameters from high-dimensional 



24 



data. Considering high-dimensional spatial observations reduces the effect 



25 



of deep uncertainty associated with different priors and results in sharper 



26 



projections based on our climate model. We demonstrate that our com- 



27 



putationally efficient approach may be widely applicable to a variety of 



28 



high- dimensional computer model calibration problems. 



29 1 Introduction 

30 Computer models play an important role in understanding complex physical pro- 

31 cesses in modern science and engineering. They are particularly important in cli- 

32 mate science where climate models, complex deterministic systems used to model 

33 climate processes, are used both to study climate phenomena as well as make pro- 

34 jections about future climate. A major source of uncertainty in climate projections 

35 is due to uncertainties about model parameters. Calibration of the parameters us- 

36 ing observational data is hence one avenue to reduce the uncertainty in future 

37 projections. A number of issues and challenges arise when performing statisti- 

38 cal calibration of model parameters. Because each run of the computer model 

39 is computationally expensive, computer model output is typically obtained for a 

40 relatively small sample of parameter values. Furthermore, the model output at 

41 each parameter setting may be high-dimensional and in the form of spatial fields. 

42 A sound statistical approach to this problem therefore needs to simultaneously 

43 address the spatial dependence in the data and model outputs, account for various 

44 sources of uncertainty, and remain computationally efficient. 

45 Previous approaches for climate model calibration have relied on heavy data 
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aggregation, turning three-dimensional data into two or one- dimensional data (e.g. 



Sanso and Forest , 2009 Goes et al. , 2010 Olson et al. , 2012) largely due to compu- 



48 tational considerations (cf. Schmittner et al. , 2009). An important and interesting 



question is what information, if any, is lost by this data aggregation. This is a 
largely unanswered question due to the inability of existing methods to analyze 
large spatial data sets of both computer model output and observations. Through- 
out this manuscript we will use "large" to refer to data sets that comprise over 
tens of thousands of spatial observations. Here we develop a computationally effi- 
cient approach that handles large data sets. This approach gives us the freedom to 
carry out a careful study of the effects of data aggregation, for example comparing 
calibration based on unaggregated three-dimensional data with calibration based 
on aggregated two-dimensional data. Our approach also enables to investigate 
the interaction between data aggregation and data-model discrepancies and errors 
when inferring computer model parameters. Our study addresses the questions of 
(i) whether aggregation is harmful or helpful, and (ii) how observations may be 
collected and processed in order to calibrate climate models and make projections 
with them. 

We adopted a Gaussian process based approach to the calibration problem 



64 (Sacks et al. , 1989 Kennedy and O'Hagan, 2001). Gaussian processes provide flex- 



ible statistical interpolators or "emulators" of the computer model across various 
parameter settings and are therefore attractive for climate model calibration (cf. 



Sanso and Forest, 2009 Bhat et al. , 2012). Unfortunately, the likelihood evalua- 



tions involved in fitting such models can become prohibitive with high-dimensional 
spatial data due to the expensive matrix operations involved. Current approaches 
for high-dimensional computer model calibration can reduce the computational 



burden and make likelihood evaluation feasible for moderately large datasets (spa- 
tial fields observed at a few thousand locations) or datasets that are on a regular 



73 and complete grid (Higdon et al. 2008 Bhat et al. , 2012; Bayarri et al. 2007). 



However, to the best of our knowledge, no current calibration approach can over- 
come the computational challenge of dealing with large spatial data sets on a 
highly irregular grid. 

The scientific problem motivating the statistical analysis is to project the future 
state of the North Atlantic meridional overturning circulation (AMOC) in response 
to anthropogenic climate forcing. The AMOC is a large-scale ocean circulation 
that transports cold and dense water equatorward in deep North Atlantic, and 
warm and salty water poleward in the upper layers of North Atlantic. The AMOC 
might show a persistent weakening in response to anthropogenic forcing. Because 
the AMOC plays an important role in heat and carbon transport, an AMOC 
weakening is projected to have considerable impacts on climate, and, in response, 



on natural and human systems (cf. Alley et al. 2007 Keller et al. 2005, 2007). 
We use previously published ensemble runs of the University of Victoria Earth 



System Climate Model (UVic ESCM) (Weaver et al. 2001) to set up an example 



es calibration problem. Vertical ocean mixing is important in projecting the AMOC 



89 (Wunsch and Ferrari, 2004), but much of the mixing occurs on scales below that 



go of the UVic ESCM, hence mixing is "parameterized" (cf. Weaver et al. , 2001 



Schmittner et al. , 2009; Goes et al. , 2010) using a "vertical background diffusivity" 



(K bg ). The AMOC projections depend on the K bg parameter values (e.g. Goes 



et al. 2010). The value of K bg is uncertain; it therefore needs to be calibrated 



94 using observations of the climate that are informative about K bg (cf. Goes et al. 



2010 Bhat et al., 2012). 



Here, we calibrate K bg using observations of ocean potential temperature from 



World Ocean Atlas 2009 (Antonov et al. 2010 Locarnini et al. 2010). K bg affects 



the depth of the oceanic pycnocline ( Gnanadesikan 1999) and the AMOC (Bryan 
1987 Goes et al. 2010). As a consequence, models with different K bg values are 



expected to result in different ocean temperature distributions. Ocean tempera- 
tures are therefore informative about K bg . We provide a detailed description of the 
data in Section |4} Note that both the UVic ESCM output and the observational 
data are spatial data sets with more than 60,000 spatial locations. Of particular 
interest is how data aggregation affects the calibration result for K bg . Often ob- 
servations of climate and climate model outputs are 3-D spatial fields. When the 
spatial data sets are large it is common practice to aggregate them into 1-D or 2-D 
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et al.| 2010 Bhat et al. 2012 ) either to avoid computational issues or because the 
skill of the models at higher resolution may not always be trusted. 

We describe a new approach that enables calibration with large spatial data 
sets without doing data aggregation. While our methodology has been devel- 
oped in the context of high-dimensional spatial fields, it applies more generally. 
In our simulated examples, we have shown that the method can handle compli- 
cated model-observation discrepancy processes without sacrificing computational 
efficiency. Our study suggests that data aggregation can cause considerable infor- 
mation loss. Reducing the data aggregation can enable valuable new insights and 
considerably reduce parametric and projection uncertainty. 

The remainder of this paper is organized as follows. In Section 2 we describe 
our two-stage framework for climate model calibration and the associated compu- 
tational challenges. In Section 3 we propose a general model calibration approach 



121 in a reduced dimensional space that uses a combination of principal components 

122 and a basis representation to overcome the computational challenges. In Section 

123 4 we describe the data and provide implementation details and in Section 5 we 

124 discuss the results from simulated examples and real data. We conclude this paper 

125 with caveats and future directions in Section 6. 



2 Model Calibration Framework 



Our computer model calibration framework consists of two stages, (i) model emu- 



128 lation and (ii) parameter calibration (Bayarri et al. , 2007 Bhat et al. , 2012). First 



129 we construct an emulator that interpolates computer model outputs at different 



i3o parameter settings using Gaussian random fields (Sacks et al. 1989). This can 



i3i be viewed as statistical interpolation or "kriging" (Cressie, 1994) in the computer 



model parameter space. Second, we infer the computer model parameters by re- 
lating observational data to computer model output using the emulator, while 
considering observational error and allowing for systematic discrepancies between 



135 the model and observations (Kennedy and O'Hagan, 2001). Note that this two 



stage approach has some advantages over fully Bayesian methods that combine the 
two stages into a single inferential step. By 'cutting feedback' or modularization 



138 (Rougier, 2008; Bhat et al. 2012 Liu et al. , 2009), this two-stage approach ensures 



that inference in the emulation stage is not contaminated by model discrepancy 
and observational error. In addition, separating the emulation stage from calibra- 
tion provides an easier way to diagnose the accuracy of an emulator. Furthermore, 
computations are faster and parameter identifiability problems are reduced. 

Let y(s, 6) denote the computer model output at the spatial location s = 



6 



144 (longitude, latitude, depth) G5CR 3 and the model parameter setting 6 © C 

145 R 3 , where S is the spatial domain of the process and is the computer model 
we parameter space. Z(s) is the corresponding observation at the spatial location s. 

147 Since each run of the climate model is computationally expensive, we can obtain 

148 computer model outputs only for a relatively small number of design points p. 

149 We denote these design points in the parameter space by 0\, . . . , p e T> C M 3 . 

150 Let Yj E W 1 be the computer model output at each parameter setting Oi for 

151 i — 1, . . . ,p. Each computer model output Yj = (Y(si, Oi), . . . , Y(s n , 0i)) T is a 

152 spatial process observed at n different spatial locations (si, . . . , s n ). Let Y be the 

153 vector of concatenated computer model outputs such that Y = (Yi, . . . , Y P ) T . We 

154 also denote the observed spatial process at n locations by Z = (Z(si), . . . , Z(s n )) T . 

155 Note that we assume that the locations for each model output and observation data 

156 are the same. If they are different, one can interpolate either of them depending 

157 on which one has a higher resolution. Our objective is to infer the parameter by 

158 combining information from Z and Y. 

159 2.1 Two-Stage Emulation and Calibration 

wo We first outline our general framework for emulation and calibration. 

lei Model Emulation Using Gaussian Random Fields 



Following Bhat et al. (2012) we approximate the climate model using a Gaussian 



process such that 

Y~ J Y0*/j,E(y), 
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162 where ftp is a linear mean function X/3 with a np x 6 covariate matrix X and a 

163 vector of regression coefficient /?. The covariates in X are the spatial locations 

164 (e.g. latitude, longitude, and depth) and the climate parameters. The covariate 

165 matrix X contains all the spatial coordinates and the parameter settings used to 
lee define the covariance matrix H(£ ). ^ y is a vector of parameters determining the 

167 covariance matrix We fit a Gaussian random field to Y by finding the 

168 maximum likelihood estimate (MLE) of (/3, £ ), denoted by (/3, £ ). 

169 The fitted Gaussian random field defines the probability model for the com- 

170 puter model output at any location s G S and parameter setting 6 e 0. Therefore, 

171 the Gaussian process model provides a predictive distribution of computer model 



172 output at any untried value of 6 given the existing output Y (Sacks et al. 1989). 

173 We denote the resulting interpolation process by rj(6,Y) and call it an emula- 

174 tor. This emulator acts as an interpolator while at the same time providing a 

175 quantification of interpolation uncertainty. 

176 Model Calibration Using Gaussian Random Field Model 

177 Once an emulator t](6, Y) is available, we model the observational data Z, 

Z = r/(Y, 0) + S + e, (1) 

178 where e ~ N(0, cr 2 I) is independently and identically distributed observational 

179 error and 6 is a data-model discrepancy term. <5 is also modeled as a Gaussian 

180 process, thus 5 ~ N(0, £d(£ d )) with a spatial covariance matrix S^^) between 

181 the locations s 1; . . . , s n and a vector of covariance parameters £ d . The details re- 



182 garding the specification of the covariance function are provided in Section 3.2.2 
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183 This discrepancy term is crucial for parameter calibration (cf. Bayarri et al. , 2007 



184 Bhat et al. 2010). Note that this problem is ill-posed without any prior informa- 

185 tion for £ d , so an informative prior is necessary. Our inference for 0, £ d and a 2 is 



186 based on their resulting posterior distribution. 



is? 2.2 Challenges with High- Dimensional Data 

las High-dimensional datasets pose considerable computational challenges due to the 

189 expensive likelihood function calculations that involve high-dimensional matrix 

190 computations. This is an important issue because model output and observational 

191 data in many fields, particularly climate science, tend to be high-dimensional. For 

192 instance, in the calibration problem described in Section |4j the dimensionality 

193 of the model output and the observational data is n = 984 in the 2-D case and 

194 n — 61, 051 in the 3-D case. The latter example involves prohibitive computations 

195 with naive implementations (discussed and explained below). For instance, with n- 

196 dimensional climate model outputs at p different parameter settings, evaluation of 

197 the likelihood function requires 0(n 3 p 3 ) operations. Therefore numerical methods 
we such as Newton-Raphson or MCMC algorithms become infeasible. For the 3-D 

199 case, this translates to 7.585 x 10 13 flops of computations, which takes more than 

200 3 hours for a high performance modern single-core processor. 
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201 3 Model Calibration with High-Dimensional Spa- 

202 tial Data 

203 We develop a dimension reduction approach based on spatial basis functions to in- 

204 crease computational efficiency. Spatial basis functions can map high-dimensional 



205 data into a low-dimensional space (Bayarri et al. 2007) and find a representation 



206 of the probability model that results in a lower computational cost for likelihood 



207 evaluations (Higdon et al., 2008 Bhat et al. , 2012). Since there may be a trade-off 



208 between parsimony and accurate inference, it is crucial to find a set of spatial ba- 

209 sis functions that gives a computationally feasible likelihood formulation without 

210 a considerable loss of information. Below, we review drawbacks to the current 
2n approaches in the context of high- dimensional spatial data and propose a new 
212 approach to overcome these limitations. 



3.1 Current Approaches 



214 A few different approaches to climate model calibration with multivariate com- 



215 puter model outputs have been developed in recent years, including Bayarri et al. 



2ie (2007), Higdon et al. (2008), and Bhat et al. (2012). These approaches, however 



217 are not readily applicable to the 3D model output and observations we consider 

218 here due to the following reasons. First, in spite of the gains in computational 

219 efficiency likelihood evaluations remain computationally prohibitive. The compu- 



220 tational cost of a single likelihood evaluation in the emulation step in |Bhat et aL 

221 (2012) scales as 0(nJ 2 ) where J is the dimensionality of the basis matrix. In 



Higdon et al. (2008), the computational cost scales with 0{p Jy) where J y is the 
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223 number of principal components used to represent the data. For the 3D calibra- 

224 tion problem we consider, n is 61,051, J is around 20,000 and J y is around 20. 

225 Second, the transformation based on the basis matrix may not be applicable to 



226 |jtwo or three-dimensional spatial data. Using a wavelet transformation (Bayarri 



227 et al. , 2007) requires the same dyadic (a power of 2) number of data points for 

228 each spatial dimension, and the data need to be on a regular grid without missing 

229 values; irregular data and missing values are common in climate model calibration 

230 problems. In addition, conducting Bayesian inference on the joint posterior dis- 

231 tribution may pose difficulties, both computationally as well as in terms of prior 



232 specification and identifiability issues. For example, Higdon et al. (2008) requires 

233 estimating 4 x J y +1 parameters, which translates to an 81-dimensional distribution 

234 for the 3-D case in Section HI 

235 3.2 Reduced Dimensional Model Calibration 

236 Our method to overcome the aforementioned challenges relies on (i) represent- 

237 ing the spatial field using principal components, and (ii) emulating each principal 

238 component separately. Instead of using principal component basis to reduce the 



239 complexity of matrix computation as in Higdon et al. (2008), we use it to map 

240 the computer model outputs into a low-dimensional space and construct an em- 

241 ulator in that space directly. Since the principal components are uncorrelated by 

242 construction, we can build the emulator by constructing a 1-dimensional Gaus- 

243 sian process for each principal component in parallel. Fitting Gaussian random 

244 fields for each principal component requires estimation of only 5 parameters. The 

245 likelihood evaluations involve covariance matrices of size p x p. These features 
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246 allow us to construct the emulator in a computationally efficient and highly au- 

247 tomated manner. Moreover, since the principal component transformation can be 

248 applied to non-dyadic spatial data with irregular locations, it has a broader range 

249 of application than wavelet transformations. In the calibration step, we develop 

250 an approach to map the observational data into a low dimensional space. 

251 3.2.1 Computer Model Emulation 

252 The first step of this approach is to find the basis matrix for computer model 

253 output. We consider the computer model outputs as an n-dimensional dataset 

254 with p replications and find the principal component basis. Let M denote the 

255 p x n matrix storing the computer model outputs Yi, . . . , Y p such that 



M 



/yr\ 



(2) 



256 Following the standard process of finding principal components, we first preprocess 

257 the computer model outputs to make the column means of the matrix M all 

258 0's. Applying singular value decomposition, we find the scaled eigenvectors ki = 

259 V^i e i ? • • • , k p = a/\^ p where Ai > A 2 > • • • > \ p are the ordered eigenvalues 

260 and ei, . . . , e p are the corresponding eigenvectors of the covariance matrix of M, 

261 where J y <C p is the number of principal components that we decide to use in the 

262 emulator. One can choose the number of principal components by looking at the 

263 proportion of explained variation given by ^T 1 / . We define the basis matrix for 

264 computer model output by K y = (k 1; . . . , kj ). 
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For each parameter setting Oi (i = 1, . . . ,p), the first J y principal components 
Y f = ( Y a > ■ * * 5 *ijJ T are computed by 



Yf = (K^K,)- 1 ^ 



265 Let Y fi = (Yf, . . . , Y R ) T , hence each element of this matrix {Y R }ij = Y R is the 

266 jth principal component at the ith computer model parameter setting. Since the 

267 columns in K y are orthogonal, the principal components found here are uncorre- 

268 lated to each other and this leads us to a parallelized emulator construction that 



269 is similar to the wavelet transformation approach in Bayarri et al. (2007). For 

270 each jth principal component across the parameter settings (i.e. for each jth col- 

271 umn of Y R ), we construct a Gaussian random field with the squared exponential 

272 covariance function such that 



CovClg, Y R ) = K y , exp ( - £ '^^ r ) + Cy,A(k = I) (3) 



273 with partial sill K y j, nugget ( y j, and range parameters (f) y j = {4> y ^j,4' y ^ji ( t ) y,i,j) T ■ 

274 Leave-10-percent-out cross-validation experiments with 50 different randomly gen- 

275 erated scenarios indicate that the squared exponential covariance shows a better 

276 fit than other alternatives such as exponential covariance. 

We denote the collection of emulator parameters for each jth principal com- 
ponent by $, y j = (n y j, ( y j, 4> y j) T . One can construct the emulator by finding the 
MLE £ ■ for each j separately. The emulator r](6, Y R ) is the collection of predic- 
tive processes of J y principal components at 6 defined by the covariance function 
([3]) and the MLEs £ y l , . . . ,£ y> j • Note that even though we constructed the em- 
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ulator in terms of the principal components, we can make a projection Y* in the 
original space at a new parameter setting 6* by computing 



Y* = K y ri(e*,Y R ). 

277 To summarize, the emulation step uses the data Yf , . . . , Yj^ of dimension p and 

278 compute MLEs £ 1; . . . , $, Jy . Hence, the computational cost is reduced from 0(n 3 p 3 ) 

279 to 0(J y p 3 ) when compared to a naive approach. The resulting fitted model is then 

280 used for the calibration step as described in the following section. 

281 3.2.2 Computer Model Calibration 

Using r)(0,Y R ), the emulator for the principal components, we define the model 
for observational data by 

Z = K y ri(0,Y R ) + K d v + e, 



282 where K^i/ is a kernel convolution representation ( Higdon , 1998 ) of the discrepancy 

283 8. v is a vector of independent and identical Normal random variates at J d <C n 

284 locations, v ~ N(0, Ij d ). ai, . . . , aj d G S. We define the kernel basis by 

etc \ - /~u~ qyi^ ( 9(sii,S2i } aij,a 2 j) \s 3i - a 3j \ \ , 

[)&4)i3 - V K <i eX P T T ' I 4 ) 

V <Pd,l 4>d,2 J 

where Ski and a^j are the kth elements of s< and aj respectively. Kd is the precision 
parameter determining the magnitude of discrepancy, and 4>d,i,<f>d,2 > are the 
range parameters specifying the bandwidth of kernels. The geodesic distance func- 
tion measures the great circle distance between two points on the Earth's surface. 
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The function g(su, s 2i , ay, a 2 ,) is given by 



r arccos (sin(s 2 j) sin(a 2 j) + cos(sin(s 2 j) cos(a 2 j) cos \su — aij\) , 



where r is the radius of Earth (6378.388km). By following Higdon et al. (2008), 



the range parameters are pre-specified by scientific expert judgment; this reduces 
computations and identifiability issues. The kernel function in ^ yields a valid 
covariance under geodesic distance since it is strictly positive definite on a sphere 



289 (Gneiting, 2011). We assumed separability for distance along the surface and dis- 



tance along the depth. The resulting process is approximately twice differentiable 



29i (Zhu and Wu, 2010), which produces a reasonable model for discrepancy. Even 



292 though the discrepancy model implies an isotropic discrepancy process (Higdon 



2002), the resulting process is flexible enough to capture the general trend in the 



discrepancy. 

Instead of using the model ^ directly, we conduct calibration with reduced- 
dimensional data for computational efficiency. Let Z R be a reduced version of the 
original data such that 



/ 



Z R = (K T K)- 1 K T Z 



T}(0,Y R ) ^ 



+ (K r K) _1 K r e, 



V " ) 



295 where K = (K y K^). The probability model of Z R is 



Z R ~ iV 



// \ / 

Mr, 



V 



V / 



Sr, 



\ \ 

+ a 2 (K T K)- 1 



\ K dlj d ) 



(5) 



/ 
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296 where fj, v and T, v are the mean and covariance given by the emulator t](6,Y R ). 

297 It is often helpful to apply singular value decomposition to and use the first 

298 c <C J d eigenvectors K^ c in place of K d to find Z R . In addition to the obvious 

299 computational advantage, this often results in better inference since it corresponds 



300 to a regularized estimate given by ridge regression (see Hastie et al. , 2009 , pg. 66); 

301 this was corroborated by our extensive simulation studies. 

302 Note that the term er 2 (K T K) _1 in ^ automatically adjusts the contribution 

303 of each principal component to the calibration result. This can be illustrated by 

304 considering the model without the discrepancy, and the variance in the likelihood 

305 function is simply T, v + cr 2 (K^K J ,)~ 1 . Since (K^K^,) -1 is a diagonal matrix and its 

306 ith diagonal element is the reciprocal of the ith eigenvalue, (K^ K y ) _1 inflates the 

307 variance of principal components with small eigenvalues. Therefore, the principal 

308 components with smaller explained variance will have less effect on the calibration 

309 result. 

310 We now briefly examine the covariance structure implied by our model. Using 
3n the leading J principal components, the covariance between computer model out- 

312 puts at two different spatial and parametric coordinates (si,0i) and (s 2 ,0 2 ) can 

313 be written as 



Cav(y(si,0i),y(s 2> 02)) ~ Gov ^^(sO^^i), ^e,(s 2 K(0 2 ) 

V i=l j=l 

J 



^2 \/^ e i( s i)v / ^i e i( s 2)Cov(w,(0 1 ),w; i (02)), 



i=i 



16 



where e»(-) is the ith eigenfunction satisfying 

Cov(y(0i,s 1 ),y(0 2 ,s 2 ))e i (s 2 )ds 1 = A^s^CovK^a), w,(0 2 )), 



with the corresponding eigenvalue A, and the process of principal component Wi(-). 
The leading eigenfunctions give the best approximation among all possible orthog- 



onal bases since it minimizes the total mean square error (Jordan, 1961). Since 
we can assume different covariance functions for each principal component pro- 
cess, our model can yield a non-separable space-parameter covariance function. In 
contrast, if we were to assume separability such that 

Cov(y(0a, si), Y(0 2 , s 2 )) = C s (s h s 2 )O,(0 1; 2 ), 

3u for some positive definite covariance functions C s and Cg, the covariance function 
315 for the ith principal component process becomes 



CovK^),^)) = Cov(y y(0 1 ,s 1 )e i (s 1 )ds 1 ,y y(0 2 ,s 2 ) ei (s 2 )ds 2 ) 

ei(si) ( [ Cov(Y(0 1 ,s 1 ),Y(0 2 ,s 2 ))e i (s 2 )ds 2 ) tfai 



C g (0i,O 2 ) J ei(si) / Cov(si,s 2 )ej(s 2 )ds 2 dsi 
Cq{0\, 2 )\i. 



316 The separability assumption therefore results in a restrictive covariance structure 

317 such that the correlation function for all principal component processes are the 

318 same. Hence even though our reduced dimensional approach utilizes a covariance 

319 that is easy to specify, it provides a richer class of covariance functions than a sepa- 
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320 rable covariance structure. Our cross-validation studies show that our assumption 

321 is adequate for emulating the computer model. 



322 Priors 

323 We estimate the joint density of 6, Ka, and a 2 using the Metropolis-Hastings al- 



324 gorithm. Following Bayarri et al. (2007), we allow for additional flexibility by 



325 estimating the partial sill parameters k V} i, . . . , K y ,j y for the emulator. Prior spec- 

326 ification for the parameters in the observational model is straightforward. The 

327 discrepancy variance k<i and the observational error variance a 2 receive inverse- 

328 gamma priors with small shape parameter values. The prior for each parameter 

329 being calibrated is a uniform distribution over a broad range or determined by sci- 

330 entitle knowledge. In order to stabilize the inference, we put an informative prior 

331 to encourage K y< i, . . . , n y ,j y to vary around their estimated values in the emulation 

332 stage. 

333 Computing 

334 The computational costs may be summarized as follows: 



335 1. Find the basis matrix K y = (\/Aiei, . . . , \/Aj^ej y ) by computing the singu- 

336 lar value decomposition of M in (|2]). This computation is of order 0(n 3 ), 

337 but it needs to be done only once. 

338 2. Compute Yr where its ith row is the transpose of (K^K y ) _1 K^Yj. 

339 3. Construct a Gaussian random field for each column of Y# by finding the 

340 MLE £ yi for each i = 1, . . . , J y . The computational cost is of order 0(p 3 ) 

341 for each likelihood evaluation. 
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342 4. Compute Z# = (K T K) 1 K T Z. The computational complexity of this step 
is 0((J y + Jd) 3 )- 

344 5. Using Metropolis-Hastings, draw an MCMC sample of 6, a 2 , and k V} i, . . . , K y ,j y 

345 from the joint posterior distribution based on the model in (|5]). The compu- 

346 tational cost for each likelihood evaluation is of order 0((J y + Jd) 3 )- 

347 The overall cost of our implementation is 0(pJy) for the emulation step and 

348 0((J y + Jd) 3 ) for the calibration step. 

349 4 Implementation Details 

350 Our goal is to build an emulator based on spatial output from UVic ESCM and to 

351 calibrate vertical ocean diffusivity (Kf, g ) using ocean potential temperature data. 

352 The UVic ESCM runs are 3-dimensional patterns of the mean ocean potential 

353 temperature over 1955-2006 at p = 250 parameter settings. The parameters con- 

354 trolling model outputs are vertical ocean diffusivity (K^g), anthropogenic aerosol 

355 scaling factor (A sd ), and climate sensitivity (C s ). Note that we converted long- 

356 wave radiation feedback factor, which is one of the original input parameters for 



357 UVic, into C s using a simple spline fit. We refer to Sriver et al. for the design 

358 points and more details of this ensemble runs. 

359 To avoid the problem related to model artifacts and sparse sampling, we ex- 



360 fluded data beyond 60°N and 80° S and 3000m in depth ( |Key et al.| |2004} |Schmit 



361 |tner et al.\ |2009| |Bhat et al.| |2012fl . UVic ESCM outputs are on a 77 (latitude) x 

362 100 (longitude) x 13 (depth) grid, but the number of locations that have non- 
363 missing observations is 65,595. The missing values occur because there is no ocean 
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at the locations in the UVic ESCM representation. 

The observational data is on a 180 (latitude) x 360 (longitude) x 33 (depth) 
grid, and we remap this observed data into the Uvic model grid using a linear in- 
terpolation This results in a relatively small reduction of data points to 61,051 data 
points. The model output locations considered in our modeling are also adjusted 
accordingly. We convert the observed in-situ temperature field into potential tem- 
perature field in order to (i) have the same measurement unit with UVic ESCM 
output and (ii) adjust the effect of pressure on ocean temperature. We obtain 



372 potential temperature from the in-situ temperature (Locarnini et al. 2010) and 



373 salinity fields (Antonov et al. , 2010) using UNESCO equation of state (UNESCO 



1981) following Bryden (1973) and Fofonoff (1977). During the conversion proce- 



375 dure, we assume a simplified ocean pressure field varying as a function of latitude 



and depth (Lovett, 1978). We apply our method to data at three different aggre- 



377 gation levels. In the 1-D case, we compute the vertical means at n = 13 different 



depth points (Goes et al. 2010). In the 2-D case, the longitudinal means are 



379 computed at n = 984 latitude and depth points (Bhat et al. 2012). We use the 



original pattern without any aggregation in the 3-D case (n = 61, 051). The num- 
ber of principal components is determined to have more than 90% of the explained 
variation. The number of components is 5 for the 1-D, 10 for the 2-D, and 20 for 
the 3-D case. We also tried to use 10 principal components for the 1-D, 20 for 
the 2-D, and 30 for the 3-D case to have more than 95% explained variation, but 
could not find any improvement in the calibration result. 

We use all 250 design points in the parameter space to build the emulator. 
We fix C s and A sc i at the default values of the UVic ESCM in the calibration 
stage and made an inference only for K^ g . The default values are 1 for A sc i and 
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389 3.81879 for C s . One may choose to integrate out these two parameters, but since 

390 the ocean temperature field lacks strong information about A sc i and C s , their es- 

391 timated posterior densities are overly dispersed. This introduces unnecessary bias 

392 in the estimate of K^g due to the highly non-linear relationship between climate 



393 parameters (Olson et al. 2012), thus we decided not to integrate out those two 

394 parameters. 



395 By following Bhat et al. (2012), we assume a flat prior with a broad range for 

396 Kf,g, from 0.05 to 0.55. The variance for the observational error (a 2 ) and the model 

397 discrepancy (k^) receive inverse-Gamma priors, and we denote them by lG(a u ,b u ) 

398 and lG(a z ,b z ). We set the shape parameters for them by a v = 2 and a z = 2. To 

399 check the sensitivity of our approach to prior specifications, we tried four different 

400 combinations, (2,2), (2,100), (100,2), and (100,100) for b v and b s . The emulator 

401 variances k Vj x, . . . , K y ,j v also receive inverse-Gamma priors with a shape parameter 

402 of 5. The scale parameters are determined to have modes at their estimated values 

403 in the emulation stage. Because of parameter identifiability problems which in turn 

404 affected computation, we fixed the range parameter for depth at 3000m and for 

405 surface as 2500km. We found that a wide range of the different settings for these 

406 parameters gave the same inference result for Kb g and hence our particular choices 

407 did not affect the results. The number of knot locations for the discrepancy kernel 

408 is 800 in the 3D case, 80 in the 2D case, and 13 in the ID case. The number of 

409 principal components used for the discrepancy is 200 in the 3D case, 20 in the 2D 

410 case, and 5 in the ID case. The number of principal components was determined 

411 using standard practice - by ensuring that at least 95% of the variability in the 

412 data was explained in each case. In order to check the robustness of our results, we 

413 tried different numbers of principal components. For example, when we increased 
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414 the number of principal components to 300 in the 3D case, 30 in the 2D case and 

415 8 in the ID case we found that we obtained virtually identical calibration results. 



416 5 Results 

417 5.1 Computational Benefit 

418 The biggest challenge in the considered analysis is the computational cost of eval- 

419 uating the likelihood function in the 3-D case, which requires dealing with 61,051- 

420 dimensional spatial datasets. To our knowledge, current approaches cannot ad- 

421 dress this problem with reasonable computational effort (discussed below). The 

422 required number of principal components is about 20 for the 3-D case for reason- 

423 able accuracy and this means we still need to invert a 5, 000 x 5, 000 covariance 

424 matrix for the likelihood evaluation in the emulation stage for the method due 

425 tO 



Higdon et al. (2008). Moreover, the number of parameters to be estimated is 



426 81, and this requires updating more than ten parameter blocks for each iteration 



of MCMC to obtain a reasonably well-mixed chain. The method of Bhat et al. 



428 (2012) requires multiplication of a J x 61, 501 matrix into another 61, 501 x J ma- 

429 trix in the likelihood evaluation with certain number of knots J, and this makes 

430 the likelihood evaluation computationally prohibitive. 

431 The computational time in emulation stage for the 1-dimensional, 2-dimensional 

432 and 3-dimensional cases for different methods are illustrated in Figure [TJ The com- 

433 puting time shown in the figure is for a system with Intel Xeon E5450 Quad-Core 

434 3.0 GHz and based on the computational complexity for matrix operations in- 

435 volved in each method. We assume that each method has been implemented using 
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436 Markov chain Monte Carlo maximum likelihood estimation with 25,000 iterations. 

437 In order to consider the dimensionality of parameter space in our comparison, we 

438 assume that five parameters are updated as a block in each MCMC run. Our 

439 method can deal with the 3-D case with a reasonable computational cost. We 

440 note that while we did not include this in our comparisons, parallel computing 

441 may further speed up computations for our approach. Since there are 81 param- 



442 eters in the approach due to Higdon et al. (2008), the method requires updating 

443 17 parameter blocks for each MCMC iteration in the 3D case. 



444 5.2 The Effect of Data Aggregation on Climate Model Cal- 

445 ibration 

446 In order to study the effect of data aggregation on climate model calibration, 

447 we conducted a study with pseudo-observational data. The simulated data are 

448 generated as follows: 



449 1. Choose the 3-D pattern of UVic ESCM output with K bg = 0.2, A sc i = 1.5, 

450 and C s = 3.975975 as the synthetic truth. 

451 2. Compute three different 3-D patterns of residuals between the observational 

452 data and the UVic model outputs with K bg = 0.1, K bg = 0.2, and K bg = 0.3. 

453 The values for A sc i and C s are the same as in Step 1. Average them over 

454 each location to compute a pseudo-residual. This is a more realistic and 

455 challenging residual than one obtained by simulation from a simple error 

456 model, for example a realization from a Gaussian process model. For brevity, 

457 we describe here just this most challenging case; our methods worked even 
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458 better in terms of posterior variance when error processes were assumed to 

459 be simpler. 

460 3. Superimpose the pseudo residual on the synthetic truth to construct pseudo 

461 observational data in 3-D. 

462 4. Aggregate the 3-D pseudo observational data into 2-D and 1-D by integrating 

463 the ocean temperature with respect to water volume. 

464 The calibration results based on this simulated example is shown in Figure [2] The 

465 sensitivity test indicates that the posterior distribution of Kb g in the 1-D and the 

466 2-D cases rely on the specification of priors. This deep uncertainty is drastically 

467 reduced when the full data set (3-D) is used. A comparison result based on the real 

468 data from Ocean Atlas 2009 is shown in Figure |3j As in the simulated example, the 

469 calibration results based on the 3-D data are more robust to the prior specification. 

470 Using the full pattern of the 3-D data has important benefits as it drastically 

471 reduces the deep uncertainty due to different prior specifications. We hypothesize 

472 that this is because the full non-aggregated spatial patterns contain more infor- 

473 mation about both the observational error and the discrepancy. In order to reflect 

474 the uncertainty due to prior choice to our density estimate for Kf, g , we show in 

475 Figure [4] the posterior distributions when the prior is assumed to be with equal 

476 probability any one of the 4 priors considered, along with the resulting AMOC 

477 projections. We define AMOC projection as the annual maximum value of the 

478 meridional overturning streamfunction in the Atlantic between 0° and 70°N. The 

479 corresponding projections for AMOC change between 1970 to 1999 mean and 2070 

480 to 2099 mean indicate that the unaggregated pattern gives a much narrower 95% 

481 predictive interval than the aggregated ones. Therefore, data aggregation increases 
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482 the deep uncertainty surrounding AMOC projection, and using unaggregated data 

483 reduces uncertainty regarding the future behavior of the AMOC. 

484 6 Discussion 

485 6.1 Summary 

486 We propose a computer model calibration approach that is computationally effi- 

487 cient even when dealing with high-dimensional data. By exploiting the orthog- 

488 onality of a principal components decomposition of the data, this method can 

489 keep the computational cost affordable for high-dimensional data with more than 

490 60,000 spatial locations and 250 parameter settings. We apply this method to 

491 the problem of calibrating a climate model parameter in an Earth System Model 

492 of Intermediate complexity based on observations of the potential temperature. 

493 The calibration results show that using 3-D data reduces the uncertainty about 

494 Kb g and is more robust to various prior specifications than calibration based on 

495 2-D or 1-D aggregated versions of the data. The results suggest that using unag- 

496 gregated data is valuable for reducing deep uncertainty associated with different 

497 priors. Our simulated examples show that our approach can handle complicated 

498 model-observation discrepancies. The method can be easily extended to allow for 

499 calibration with multiple tracers - we can simply consider the variance-covariance 

500 matrix for all tracers and use its principal components to build an emulator. 
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501 6.2 Caveats 

502 A general issue with principal components is also worth considering in this con- 

503 text: the principal components for the computer model outputs are selected based 

504 on explained variation, and thus there is no guarantee that these leading principal 

505 components carry the most important information about the climate parameters. 

506 However, our extensive study of the effect of changing the number of principal 

507 components suggests that this is not problematic in our context. Our results are 



508 consistent with the recent theoretical results in Artemiou and Li (2009) that sug- 

509 gest that there is a low probability that other (non-leading) principal components 

510 will have a strong correlation with the climate parameters. We hypothesize that 

511 our principal components based approach does not lose valuable information about 

512 the climate parameters. In the discrepancy model, one important simplifying as- 

513 sumption is the separability between surface and depth effects. Our simulated 
5u example shows that the separability assumption provides a good approximation 

515 to the realistic discrepancy process. Non-separable covariance function that com- 

516 bines geodesic distance and Euclidean distance remains as an avenue of ongoing 

517 and alive research and subject of future work. Furthermore, our study of calibra- 

518 tion with simulated examples shows that even though the number of Kf, g settings at 

519 which the model is run is relatively sparse, there is enough information to reliably 

520 calibrate K^g based on our emulator. 

521 Our study is also subject to usual caveats with respect to scientific conclu- 

522 sions. First, we ignore the interpolation uncertainty when we compute the density 

523 of AMOC projection based on the density of Kj, g . Second, the result is based on 

524 a single data set, and thus we cannot fully evaluate the effect of structural uncer- 
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525 tainty due to the model-observation discrepancy and unresolved natural variability 



526 ^annot be accounted for; this variability could impact conclusions as well (Olson 



et al. ). These caveats, of course, apply to almost all existing approaches to climate 



528 model calibration and projection. 
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Figure 1: Comparison of computational costs for the emulation step between the 
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Figure 2: Prior sensitivity test in the simulated example. Calibration of K^ g 
value based on (a) 1-D depth profile (upper panel), (b) 2-D latitude-depth pattern 
(middle panel) (c) 3-D non-aggregated data (lower panel). Each line represents 
posterior density from four different priors: (b v = 2, b z = 2) (Prior 1, solid line), 
(b v = 2 } b z = 100) (Prior 2, dashed line), (b v = 100, b z = 2) (Prior 3, dotted 
line), and (b v = 100, b z = 100) (Prior 4, dot-dashed line). The solid vertical line 
represents the true value of Kj, g in the synthetic truth. 
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Figure 3: Prior sensitivity test using observational data from the World Ocean 
Atlas 2009. Calibration results based on (a) 1-D depth profile (upper panel), (b) 2- 
D latitude-depth pattern (middle panel) (c) 3-D non-aggregated data (lower panel). 
Each line represents posterior density from four different priors: (b v = 2, b z = 2) 
(Prior 1, solid line), (b v = 2,b z = 100) (Prior 2, dashed line), (b v = W0,b z = 2) 
(Prior 3, dotted line), and (b v = 100, b z = 100) (Prior 4, dot-dashed line). 
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Figure 4: Combined posterior densities of Kf, g from different prior specifications 
(lower left), the relationship between Kb g and projected AMOC change of the 2070 
- 2099 mean from the 1970 - 1999 mean (upper left), and the resulting AMOC 
change projections (upper right) using 1-D (solid line), 2-D (dashed line) and 3-D 
(dotted line) data with their 95% credible intervals (bars at the right) 
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